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ABSTRACT 

We simulate the growth of isolated dark matter haloes from self-similar and spherically sym- 
metric initial conditions. Our N-body code integrates the geodesic deviation equation in order 
to track the streams and caustics associated with individual simulation particles. The radial 
orbit instability causes our haloes to develop major-to-minor axis ratios approaching 10 to 1 
in their inner regions. They grow similarly in time and have similar density profiles to the 
spherical similarity solution, but their detailed structure is very different. The higher dimen- 
sionality of the orbits causes their stream and caustic densities to drop much more rapidly than 
in the similarity solution. This results in a corresponding increase in the number of streams 
at each point. At 1 % of the turnaround radius (corresponding roughly to the Sun's position in 
the Milky Way) we find of order 10^ streams in our simulations, as compared to 10^ in the 
similarity solution. The number of caustics in the inner halo increases by a factor of several, 
because a typical orbit has six turning points rather than one, but caustic densities drop by 
a much larger factor This reduces the caustic contribution to the annihilation radiation. For 
the region between 1 % and 50% of the turnaround radius, this is 4% of the total in our sim- 
ulated haloes, as compared to 6.5% in the similarity solution. Caustics contribute much less 
at smaller radii. These numbers assume a 100 GeV c^^ neutralino with present-day velocity 
dispersion 0.03 cm s^^, but reducing the dispersion by ten orders of magnitude only doubles 
the caustic luminosity. We conclude that caustics will be unobservable in the inner parts of 
haloes. Only the outermost caustic might potentially be detectable. 
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1 INTRODUCTION 

Dark matter is thought to consist of weakly interacting particles 
that are cold, meaning that their thermal velocities in unclustered 
regions of the present Universe are very small. For example, for 
standard values of the relevant cross-sections, a neutralino with 
a mass of 100 GeV is predicted to have a velocity disper- 
sion of just 0.03 cm/s in such regions, corresponding to a very 
high phase-space density. Since the particles are collisionless, this 
phase-space density is conserved along particle trajectories, and 
the nonlinear collapse of dark haloes gives rise to caustics. For in- 
finitely cold dark matter, the density is formally infinite at caus- 
tics, because the particles occupy a three-dimensional "sheet" in 
six-dimensional phase-space and projection of this sheet onto con- 
figuration space leads to singularities. These singularities are regu- 
larised in realistic cases by the small but finite velocity dispersion 
of the dark matter particles. These caustics were first discussed by 
Arnold et al. (1982) and ZePdovich et al. (1983) for the case of 
a neutrino-dominated universe. Their properties were worked out 
in simple, one-dimensional similarity solutions for the growth of 
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sph erical haloes by Fillmore & Goldreicri im4) (FG hereafter) 
and lBertschingeJ jl985t) . 

It has been suggested that caustics might have an impact 
on various observables, in particular gravitational tens ing (e.g. 
Hogan 1999; Charmousis et al. 2003; Gavazzi et al. 2006) or anni- 
hilation radiation (e.g. Hogan 2001; Mohavaee & Shandarinl 



2006; Moha vaee et al J l200l iNataraian & Sikivid |2008| ; 
Moh avaee & Salatil |2008| ; lAfshordi et alj |2009|) . Over the last 
10 years most studies of the properties of caustics have been based 
on spheric al models fo r self-similar infall such as those of FG 
and ,Bertschingeil ( Il985l) . or extensions where a ngular momentum 
is in tr oduced to avoid pu r ely radial orb it s jWhite & Zaritskvl 
'1992"; 'Sikivie et al.' '1995"; 'Nusser' 1200 ll ; INataraian & Sikiviel 
2006; Natarajan 2007; Sikivie 1998). As we will see below, the 
original similarity solutio ns are violent ly unstable to the radial 
orbit instability (ROI) tAntonovl Il973h . Exact density profiles 
for dark matter caustics w e re firs t calculated for the FG model 
by Mohavaee & Sha ndarinl ( l2006h (MS hereafter) following an 
approach originall y developed by jZeldovich & Shan d arin ( 198^) 
and iKotok & Shandariij ( 1 19871) . In contrast, iHoganI feoOlh dis- 
cussed the effect of caustics on annihilation radiation using general 
arguments which avoid simplified halo models, estimating that 
caustics might boost the annihilation luminosity by a factor ~ 5 in 
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the outer halo, a significantly larger factor than found by MS for 
the spherical similarity solution. 

Although very large simulati ons of the formation of dark 
haloe s are now feasible (e.g. Die mand et alj|200^ : ISpringel et al.l 
l2008l;IZemDetai]|2009h . it is still not possible to resolve caustics 
using standard N-body techniques. Here we us e the geodesic devi- 
ation equation (GDE) formalism presented in IVogelsberger et al.l 
12008) (VWHS hereafter) to gether with a rigorous treatment of 
caustic passages derived in IWhite & Vogelsbe rger (2009) (WV 
hereafter) to follow the evolution of caustics in fully three- 
dimensional simulations of halo formation. For simplicity, we start 
from self-similar, spherical initial conditions, but the ROI causes 
our simulated haloes to develop into highly elongated bars with a 
detailed structure quite different from that of the spherical simi- 
larity solution. This has substantial implications for the number of 
dark matter streams predicted near the Sun and for the importance 
of annihilation radiation from caustics. 

The plan of our paper is as follows. In Section 2 we de- 
scribe our initial conditions and the numerical techniques used 
for our simulations. In Section 3 we demonstrate that our N-body 
implementation is working correctly by solving a simple spheri- 
cal test problem. In Section 4 we discuss results from our three- 
dimensional simulations, contrasting them with the predictions of 
the similarity solution. We begin with the shape and the density 
and velocity dispersion profiles, and we move on to fine-grained 
phase-space structure and the caustic annihilation rate. We give our 
conclusions in Section 5. 



2 NUMERICAL TECHNIQUES 
2.1 Initial conditions 

We start our simulations from self-similar initial conditions, where 
most of the particles are in the linear regime, i.e. well beyond 
the initial turnaround radius. The similarity solutions assume 
an Einstein-de Sitter universe in which the linear mass pertur- 
bation 5 Mi within a sphere containing unperturbed mass Mi, 
when extrapolated to the initial time tinitiai, satisfies 5Mi/Mi = 
1.0624(Mi/Mo)~'^, where e is a scaling index and A/q is a refer- 
ence mass taken to be the mass within the turnaround radius at the 
initial time. The equations of motion for a particle with radial dis- 
tance r and radial velocity v can then be cast into similarity form 
by introducing the variables A = r/rta and r = t/tta, where 
Ml'""*' and tta oc M^'^^ are the turnaround radius and 
turnaround time of the particle under consideration. The solution 
A(r) of the resulting equation then fully describes the phase-space 
structure of the halo at all times (see the Appendix and FG for more 
details on these similarity solutions). Note that the turnaround ra- 
dius rta can be viewed either as a function of enclosed mass, as 
here, or as a function of time, as will often be the case in the fol- 
lowing. 

We set up our initial conditions in the following way. As a first 
step we create a uniform gravitational glass within a cubic box with 
periodic boundary conditions (White 1996). We then cut out the 
largest sphere contained within the box and use it to represent an 
unperturbed Einstein-de Sitter universfl We construct similarity 



^ Note that this cut means that a cubic glass with A'^'^ particles will pro- 
duce a sphere with N'^ /(6/n) = pailicles. In the following we will 
always specify the resolution of our simulations by giving the number of 
particles in the original cube. 
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Figure 1. Top panel: Time evolution of the caustic radii in the softened sim- 
ilai'ity solution (e = 2/3, rj = 0.15) simulated with a (one-dimensional) 
shellcode. Different colours represent different caustics from red (outer- 
most) to grey (innermost). We note that the softened similarity equation 
only produces a finite number of caustics while the original FG equations 
produce an infinite number when approaching the halo centre. The sphere 
radii grow linearly in time because the turnaround radius increases linearly 
in time for e = 2/3. Bottom panel: Slice through the caustic spheres. The 
thickness of the slice is 0.0025 rta(t)- We overlaid all caustics for all times 
using the similaiity scaling. This produces exact spheres (rings in this slice). 
These plots demonstrate that caustics ai'e tracked correctly over the full sim- 
ulation period by our GDE method. 



initial conditions at time tinitiai by mapping this uniform distribu- 
tion r to the desired initial state by scaling the radial coordinate 
and setting velocities according to the similarity solution A(r) 



r(tinitial) = rta(Mo 
"(tinitial) = 



A(r) 
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na(AJo ) dA/dr 

iinitial 



^glass 



-l/3+2/(9e) 



'■gla 



(1) 



where the similarity time variable r is determined from the en- 
closed mass according to r = {AIi/Mo)~^'^^^. 

We run our simulations from tinitiai to to ~ 2/(3 -ffo) — 
9.1 Gyr, where we adopt Ho = 72 km Mpc^^ as the Hubble 
constant today. If we require that the turnaround of the bounding 
sphere occurs at to, choosing the scaling index e and the mass frac- 
tion Mo/Mtot within the initial turnaround radius determines the 
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Figure 2. Maximum caustic densities in units of tlie stream density at 
turnaround po ■ Black dashed lines show maximum densities for the analytic 
solution. The densities increase with time, because the velocity dispersion 
tTi, decreases with time. 

initial time as tinitiai = to{Mo/AItotf"^^'^ and defines a unique 
mapping from the uniform spherical glass to the similarity solution 
at tinitiai- In the following we will always use this kind of initial 
conditions. 



2.2 Simulation code 

For our simulations we adapt the GAD G ET- 3 code that was d evel- 
oped originally for the Aquarius project jSpringel et al.ll2008h . We 
apply vacuum boundary conditions, use (unless otherwise noted) a 
spline softening kernel with constant smoothing length in comov- 
ing coordinates (oc t^^^ in an Einstein-de Sitter universe), and run 
our simulations in physical coordinates. We modified the code to 
integrate the GDE as described in VWHS and to track radiation 
from each particle due to annihilations within its own fine-grained 
phase-space stream, as described in WV. This automatically ac- 
counts for the enhanced radiation as particles pass through caus- 
tics. We set the GADGET force accuracy para meter to 10~^ and 
chose a time integration accuracy of 10~^ (see lSDringelll200j . for 
the significance of these terms). We checked that these settings are 
sufficient to produce reliable results. 

To gain performance and accuracy we integrate the GDE di- 
rectly for each particle only after it has passed through turnaround. 
At turnaround all variables related to the fine-grained phase-space 
density are set to the values predicted by the analytic FG solution. 
We implement this by checking in each drift operation whether 
a given particle is currently turning around and initialising the 
integration of the phase-space distortion tensor when this is the 
case. We demonstrate below that using the analytic solution until 
turnaround is well-justified. As a particle turns around at tta, its 
phase-space distortion tensor is set to unity and its stream den- 
sity, i.e. the local 3-density of the particular stream in which it 
lives, is set to the dark matter density of the similarity solution 
at the turnaround point po = 97r^/(16(3e -I- 1)) pb(ita), where 
Pb(ita) = l/(67rGtta) is the mean density of the Einstein-de Sit- 
ter background at time fta- The (fine-grained) dark matter velocity 
dispersion ctq at that time and position can be calculated using con- 
servation of phase-space density po/o"o = Pb(io)/o"b(to)"^, where 



(Tb(to) is the dark matter velocity dispersion predicted for unclus- 
tered matter today. For a neutralino of mass nip we have (Jb (to) ~ 
lO^^^c (GeV c~^/mp)^''^. Unless otherwise noted, we will as- 
sume a 100 GeV neutralino, resulting in (7b(to) = 0.03 cm/s. 
Finally, the sheet orientation at turnaround is set to 
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4 J 3 + 1/etta' 

where Vi (q) are the three components of the initial cold dark matter 
sheet (see VWHS and WV). All GDE simulations presented in this 
paper have the fine-grained phase-space quantities initialised for 
each particle at turnaround in this way. 



3 RESULTS 

3.1 Embedded ID tests 

To test the GDE implementation in GADGET- 3 it is desirable to 
have a system where we can analytically evaluate the fine-grained 
phase-space structure. This is straightforward for the FG similar- 
ity solution (see the Appendix for details). Starting with the ini- 
tial conditions described above we might expect to be able to re- 
produce the similarity solution with our code. This is not pos- 
sible, however, because the FG solution is violently unstable to 
both radial and non-radial perturbations, as we show below. Al- 
though this result does not appear to have been demonstrated be- 
fore for the specific case of the FG solution, it is expected given 
earlier work on other spherical collapse models. This has found in- 
stabilities not only in the fully 3-dimensional case where the ra- 
dial orbit insta bility (ROI) tur ns the system into a highly elon- 



gated bar (e.g. Antonovl 19731: IP olvachenko 1 98ll: Isamej 



Carpintero & Muzzid ll995l : iMerritt., 199ft : ,MacMillan et alj 



198^ : 



20061 : 



Bellovarv et al. 2008), but also in ID (i.e. enforcing spherical sym- 
metry and radial orbits) where the regular phase-space p attern of 
the s i milari ty solution still gets destroyed tHenriksen & Widrowl 
Il99lll999h . These instabilities have a dramatic impact on the fine- 
grained phase-space structure as we will show below, so we cannot 
follow this route to check our standard code against the known sim- 
ilarity solution. 

A test can nevertheless be carried out by keeping the calcu- 
lation artificially stable. It turns out that the following approach 
works. We replaced the GAD GET- 3 treecode by a shellcode, 
where radial forces are calculated based purely on the enclosed 
mass. Particles in the initial conditions then represent mass shells. 
This removes the degree of freedom exploited by the ROI, but does 
not stabilise the system against purely 1-D instabilities. To avoid 
the destruction of the similarity phase-space pattern, we must in 
addition soften quite strongly the gravitational potential. One can 
show that a Plummer-softening that scales with the turnaround ra- 
dius still allows a similarity solution. The similarity equation is then 
a slight extension to the FG equation 



dr2 



X 



(A2 + (^A)2; 

where A4 is the dimensionless enclosed mass 



M 



2 

3^ 



^1+2/(36) 



H 



A(e) 
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H the Heaviside function, A(r) = r^/^+^Z^^"^ and 77 is the 
Plummer-softening length in units of the turnaround radius. We 
can solve this softened similarity equation in the same way as the 
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Figure 3. Top panel: Phase-space portrait at to for the softened similarity 
system. The various colours indicate how many caustics each particle has 
passed. Infalling pailicles have seen no caustic (black line). They go through 
the centre, and move out to the first apocentre, where they pass the outer- 
most caustic. The colour changes to red at this point and stays red until the 
second caustic is reached, and so on. We note that the softened similarity 
solution has only a finite number of caustics. This is why the phase-space is 
not densely occupied at low r/rta as in the case for the original unsoftened 
FG solution. Bottom panel: Densities of individual dark matter streams in 
the interval r/rta S (0.1, 0.105). At this distance the portrait shows nine 
streams (the innermost stream is coloured in yellow and particles belong- 
ing to it have passed four caustics). In black we show a histogram of their 
densities in the GDE simulation and in red the analytically predicted stream 
densities. The agreement is good, demonstrating that our method can cor- 
rectly recover the density of individual streams. 



original FG equation, including the parallel GDE integration as de- 
scribed in the appendix. With this setup we can check our N-body 
implementation by comparing the results to the analytic solution of 
the softened similarity equation. 

The following results are based on a TV = 128^ simulation 
with this method. We have chosen e = 2/3, rta(to) = 1400 kpc, 
^initial ~ ^o/lO^ and Tj — 0.15, where rta denotes the global 
turnaround radius of the halo. 

Fig. [T] shows how well our code reproduces the caustic struc- 
ture of this test problem. The top panel shows the time evolu- 
tion of the caustic radii. They grow linearly in time, because the 
turnaround radius is proportional to t for e = 2/3. The different 



colours represent the different caustic spheres (red is the outermost 
sphere, green the second, and so on). We note that the softened FG 
equation has only a finite number of caustics for large r) (nine for 
ri — 0.15) while the original FG solution (the 77 — > limit) pro- 
duces an infinite number of caustics. The bottom panel shows a 
slice through the caustic spheres from a wide range of times after 
scaling each to its own turnaround radius. In this plot we thus over- 
lay the complete time evolution of the system. One can clearly see 
that the caustics build perfect spheres that scale exactly as expected. 

The maximum densities in these caustics are shown in Fig. [2] 
The agreement of the GADGET- 3 shellcode results with the ana- 
lytic solution (black dashed lines) is good. The maximum density 
Pmax/po is proportional to 1/ ^ ab{t) for e — 2/3, and therefore 
increases with time due to the decrease of the velocity dispersion. 

In Fig. [3] we focus on the density of individual fine-grained 
dark matter streams. We determine the density of individual 
streams by binning ps,i/po,i for all particles i in r/rta G 
(0.1,0.105), where ps,i is the density of the stream particle i is 
embedded in and po,i is its stream density at turnaround. From the 
phase-space portrait Fig. [3] (top panel) one can see that there are 9 
distinct streams in this radial interval. In the bottom panel we show 
the resulting stream density histogram in black. In red we over- 
plot the analytic result for the densities of the nine streams. The 
agreement is good, showing that we correctly recover the density 
of individual streams in this system. 

We will now briefly described how we calculate the intra- 
stream annihilation rate for each individual simulation particle, that 
is the annihilation rate due to encounters with other particles in 
its own DM stream (see VWHS, WV). Each simulation particle 
(mass rrii) represents many dark matter particles (mass nip). The 
intra-stream annihilation rate of one dark matter particle is given 
by [{oAv) /'mp)pa,i. Therefore, the annihilation rate of a simula- 
tion particle is given by {{aAv) /mp){ps,iini). To get an estimate 
of the intra-stream annihilation rate, we integrate the stream den- 
sity along the trajectories of all particles as the simulation is run. 
This yields for every particle and at each snapshot time tk the time- 
integrated rate 



A^{tk) 



dt ps,i{t), 



(5) 



where we set the particle physics prefactor (aAv) /nip to unity. 
When passing through a caustic it is necessary to calculate the time 
integral analytically as described in WV, since the simulation time- 
stepping is usually not fine enough. Based on Ai{tk) and Ai{ti) 
(tk < ti) we can calculate the intra-stream annihilation rate of par- 
ticle i 



(-fi)intra — rllj 



Ajjtl) - A^{t^:) 

At 



(6) 



In the following we apply this scheme to calculate the intra-stream 
annihilation rate of individual particles. 

3.2 3D results 

In this section we study halo formation from the same initial con- 
ditions as above, but with fully 3D gravity. Thus we replace the ID 
shellcode of the last section with our modified version of the tree- 
code GADGET-3. In addition, we choose a much smaller softening. 
The ROI then acts with full force and the system develops a highly 
elongated bar in its inner regions. Our goal here is to study how 
this more complicated dynamical structure affects the fine-grained 
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Figure 4. Projected density of a slice of thickness 0.25rta through the 3D 
simulation at the final time when rta equals the radius of of the entire sim- 
ulated region. The strong bar is clearly visible and is oriented to lie in the 
plane of the slice. It formed through the action of the ROI. The red region 
is bounded by the outermost caustic. 




Figure 5. The turnaround radius evolution of the 3D simulation is compared 
to the expectation from the ID analytic similarity solution. Although the 
structure of the 3D halo is very different from that of the similarity solution, 
the turnaround radius grows in exactly the same way, linearly in time for 
e = 2/3. 

phase-space of the system, and what implications this has for the 
intra-stream annihilation rate. In the following we use a N — 128^ 
simulation with a comoving softening length of 0.01 rta(fo). We 
checked that our results do not depend on resolution (either particle 
number or softening) in any significant way. 

In Fig.|4]we show the projected density of a slice through the 
centre of the halo at to with a thickness of 0.25 rta(to)- At this 
time the last simulation particles are just reaching turnaround. The 
blue region marks the halo beyond the outermost caustic, hence the 
one-stream regime. The change to red marks the region of mul- 
tiple streams. The occupied region has a sharp spherical edge at 
rta, showing that the influence of our vacuum boundary condition 
has not propagated into the inner regions of interest. We checked 



this explicitly with other simulations, and also by verifying that the 
results given below are very similar to those at 0.5to once these 
are scaled up according to the similarity scaling. Clearly visible in 
Fig. |4] is a strong bar produced by the ROI. To show it optimally, 
we have oriented our slice perpendicular to its minor axis. At the 
time, shown we estimate axis ratios b/a = 0.41 and c/a = 0.22 
from the moment of inertia of all particles within rta, axis ratios 
b/a = 0.34 and c/a = 0.17 from all particles within 0.5 rta, and 
b/a = 0.28 and c/a — 0.12 from those within 0.1 rta- These 
numbers are quite similar att = 0.5to suggesting that the system 
is evolving in a nearly self-similar, though non-spherical fashion. 

Although the inner structure of the halo is very different from 
the FG model, the turnaround radius grows in time just as this 
model predicts. This is demonstrated in Fig. |5] where the dashed 
black line is the FG prediction for e — 2/3. This behaviour is ex- 
pected since the growth of the turnaround radius depends on the 
monopole term in the halo mass distribution and this is not affected 
by non-spherical contributions from the inner part. We note that 
this is important for our GDE integrations because we use the ana- 
lytic FG solution to initialise the GDE variables at the turnaround 
radius. 

It is obvious that the strong departures from spherical symme- 
try must nevertheless substantially change the structure of the sys- 
tem. The changes are less dramatic than might be expected, how- 
ever. In Fig. [6] we plot the spherically averaged density profile as a 
function of radius at to both for the 3D simulation and for the sim- 
ilarity solution. Except for the caustic spikes, the mean halo den- 
sity profile in agrees very well over most of the plotted radial range 
with the similarity prediction. The deviations visible on the smallest 
scales in Fig.|6]are an effect of force softening. We have performed 
simulations with 20 times smaller softening, finding that the spher- 
ically averaged profile does not deviate significantly from isother- 
mal form (p oc r^^) down to 10^'^ rta- We have also checked that 
if we change the value of the similarity exponent e in our initial 
conditions, the inner profiles of the resulting bar-like haloes are al- 
ways well described by a power law with exponent 9e/(l + 3e), 
the val ue predicted by spherical similarity solutions with nonrad ial 
orbits l lWhite & Zaritskvll992l : rSikivie et al.lll995l : lNusseJl200lh . 

The ROI disturbs the overall radial structure of the system very 
little. We see no sign of a tendency to drive the system towards a 
"unive rsal" NFW-like profile of the kind which iMacMillan et al.l 
( l2006h and lBellovarv et aT found in their own experiments 

from non-similarity initial conditions. This demonstrates that the 
ROI does not of itself produce NEW structure, though it may, of 
course, be acting in concert with other processes in these earlier 
experiments. 

The aspherical structure of the inner regions induces non- 
radial motions in the infalling material, producing a phase-space 
distribution which is six-dimensional rather than two-dimensional 
as in the similarity solution. We illustrate how this affects the ve- 
locity dispersion structure in Fig. |7] The upper panel shows total 
and radial velocity dispersions (relative to mean radial streaming) 
as functions of radius in units of rta, while the lower panel is a 
similar plot for the velocity anisotropy parameter /9(r). In both pan- 
els the 3D simulation result (which is at to) is compared with the 
spherical similarity solution. The total velocity dispersion of our 
3D halo tracks that in the similarity solution quite accurately down 
to 0.04rta, once the features due to the spherical caustics in the 
latter are smoothed over. At smaller radii it turns over and drops 
significantly. This reinforces the conclusion from Fig. [6] that the 
coarse-grained structure of our simulated halo is similar to that of 
the similarity solution despite their difference in shape. While at 
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Figure 6. Spherically averaged density profile of the 3D simulation com- 
pared to the FG similarity solution. Apart from the clear caustic spikes in 
the ID case, the two density profiles agree very well despite the very large 
shape difference between the corresponding objects. The small deviation at 
r < O.OlSrta is due to the softening of the simulation. 



r > 0.05rta orbits in our bar-like halo are predominantly radial 
and so resemble those in the FG similarity solution, at smaller radii 
the velocity distribution becomes much closer to isotropic and the 
radial velocity dispersion begins to decline towards the centre. 

The ROI leads to a considerably more complex fine-grained 
phase-space pattern in our 3D simulation than in the FG similarity 
solution. Phase-space is two-dimensional in the latter, with the par- 
ticles occupying a one-dimensional subspace which is fixed in sim- 
ilarity variables. In contrast, phase-space is six-dimensional in our 
simulation, with the particles occupying a heavily wrapped three- 
dimensional subspace. When this subspace is projected onto the 
two-dimensional phase-space of the similarity solution, the fact that 
the particles occupy a low-dimensional subspace is no longer evi- 
dent. We show this explicitly in Fig. [8] where the pattern of the 
ID similarity solution (the grey lines) is compared to the projected 
particle distribution at to in our simulation. The colours in the 3D 
case mark the number of caustics each particle has passed. As in 
the ID case, this number increases towards the centre, because of 
the shorter orbital periods at smaller radii. Near and beyond the 
turnaround radius particles follow the similarity solution, but devia- 
tions are already visible at r ~ 0.5 rta, where the infalling particles 
no longer lie exactly on the analytic sheet. 

The increase in dimensionality of the phase-space distribu- 
tion has a dramatic effect on the density of individual dark matter 
streams. Away from caustics the stream density surrounding a par- 
ticular particle is expected to decrease from its value at turnaround 
in proportion to (t/tta) ^ for an effectively one-dimensional sys- 
tem like the similarity solution, but in proportion t o (t/tta)~^ for 
a thre e-dimensional system like our simulation iHelmi & White! 
1 19991 and VWHS) . In Fig.|9]we show the median of the normed 
stream density (i.e. the stream density relative to its value at 
turnaround) for particles in radial bins at t = to, again compar- 
ing results for the simulation and for the ID similarity solution. 
Beyond the outermost caustic, stream densities are very similar in 
the two cases. At smaller radii the stream density dilution increases 
much faster towards the centre in the simulation than in the similar- 
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Figure 7. Top panel: Total (cr^ = ct^ + cr^, where Ur and at are the 
velocity dispersions in the radial and tangential directions, respectively) and 
radial (c^) velocity dispersion profiles for the 3D halo are compared to 
the radial velocity dispersion profile of the FG similarity solution. Bottom 
panel: Velocity anisotropy profile /3 = 1 — cr^/ (2(7^). Apart from caustic 
features the total velocity dispersion behaviour is very similar in the two 
models for r > 0.04rta. In the outer part of the halo, orbits are primarily 
radial (/3 ~ 1), but the velocity distribution becomes more nearly isotropic 
for r < 0.05rta and this causes both dispersions to drop in the inner regions 
as in spherical similarity solutions with nonradial orbits. 

ity solution. In the former case the inner behaviour is very close to a 
power-law. (Recall that for both models the typical orbital period at 
each radius is roughly proportional to radius, so a power law close 
to is expected in the ID case.) In the simulation the variation 
of dilution with radius is less smooth, but is still moderately well 
represented by a power-law which, as expected, is the cube of the 
one which best fits the similarity solution. 

Fig. [To] presents the radial variation of stream density in a 
slightly different way. The stream density associated with each par- 
ticle is here divided by the current mean density of the Universe, 
rather than by its value at turnaround. The plot is then constructed 
by taking medians in 50 equal logarithmic bins in radius, just as 
in Fig. |9] In the ID similarity solution, stream densities increase 
quite strongly (approximately as r~^) towards the centre, because 
the dilution effect seen in Fig.|9]is more than compensated by the 
fact that particles near the centre typically turned around earlier. 
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Figure 8. Phase-space portrait at to for the 3D simulation and (overplotted 
in grey) for the ID similarity solution for e = 2/3. The ROI destroys the 
clear phase-space pattern seen in the similarity solution. The coarse shape 
in the two cases is similar, but the 3D solution does not have a simple fine- 
grained structure in this particular projection of its 6-dimensional phase- 
space. For the 3D solution the colours encode the number of caustics a 
given particle has passed (see colourbar). One can clearly see that this grows 
towards the centre, exceeding 200 in the innermost regions. 




Figure 9. Median normed stream density (i.e. the stream density surround- 
ing each particle in units of its value at turnaround) for the ID similarity 
solution (the thick blue line) and for the fully 3D simulation (the thick red 
line). The ID relation was calculated from the exact (non-softened) simi- 
larity solution for e = 2/3. Due to the fully three-dimensional structure 
of the orbits, stream densities are diluted much faster towards the centre 
in the simulation than in the similarity solution. In both cases the median 
was calculated for 50 equal logarithmic bins in radius. This smooths out the 
caustics of the ID solution except for the outermost few. The thin blue line 
is a power-law fit ps oc r"'* to the normed stream density in the inner part 
of the similarity solution. The blue dashed line is proportional to the ID 
density dilution to the power of three and, as expected, agrees fairly well 
with the sfi'eam dilution in the inner regions of the 3D simulation. 
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Figure 10. The median stream density in units of the current cosmic mean 
density of the particles in each of 50 logarithmic bins in radius at time to . 
Results are shown for the ID similaiity solution (blue) and for the 3D sim- 
ulation (red). While the stream density increases towards the centre for the 
similarity solution, it is almost constant with radius in the simulation. 




Figure 11. Median number of caustic passages experienced by particles in 
a set of logarithmic bins in radius. The result for the ID similarity solution 
is shown as a solid blue line, that for the 3D simulation in red. As expected 
due to the more complex orbit structure, the number of caustics at a given 
radius increases by about a factor of three from ID to 3D. The bump near 
r/rta. = 0.1 is related to the complex feature in phase-space seen in Fig.[8] 

and so had higher stream densities at turnaround. For the 3D sim- 
ulation, on the other hand, the much stronger dilution towards the 
centre results in typical stream densities which are approximately 
constant with radius at any given time. Note that both in this figure 
and in the last, medians are taken over the stream density distribu- 
tion of the particles in each radial bin. Thus they give the stream 
density value that splits the mass at each radius in half. As we will 
see below, this is much larger than the density of a typical stream 
at that radius. 

Typical stream densities are quite similar near the outermost 
caustics in the ID similarity solution and in the 3D simulation, but 
they become very different in the inner regions. In contrast, the 
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Figure 12. Number of streams as a function of radius for tlie ID similar- 
ity solution (blue) and the 3D simulation (red). In 3D. stream stretching 
reduces individual stream density more efficiently than in ID, leading to 
a larger number of streams at each radius. For the ID similarity case the 
number of streams in the FG solution is shown. In the 3D case we use the 
mean harmonic stream density of the particles in a set of logarithmic bins 
to estimate the number of streams (see text). 
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Figure 14. Local ratio of the intra-stream annihilation rate to the smoothed 
annihilation rate as a function of radius. The smoothed annihilation rate is 
calculated from an SPH-based density estimate using 64 neighbours. As 
expected, the strongest caustic contributions are found in the outer regions. 
Near 1% of the turnaround radius the contiibution from caustics is at the 
percent level. This reflects the rapid dilution of stream densities in a 3D 
system. 




Figure 13. Densities at caustic passage as a function of radius and in units of 
the current mean cosmic density for the ID similarity solution (blue crosses) 
and for the 3D simulation (red curves showing the median and quartiles 
at each radius). These caustic densities are compared to the mean density 
profiles of the two systems (blue curve for the ID similarity solution, green 
curve for the 3D simulation). In 3D the outermost caustics are as dense as 
those found in ID, reflecting the more rapid mixing in the three-dimensional 
system. 



typical number of caustics varies with radius in the same way in 
the two cases. This can be seen in Fig. [TT] where we plot the me- 
dian number of caustic passages experienced by particles in a set of 
logarithmic radial bins. The number is typically three to six times 
larger for the simulation than for the similarity solution. This is eas- 
ily understood as reflecting the increased number of turning points 
along typical orbits. The overshoot in the regions just inside the 



outermost caustic is related to the feature in phase-space seen in 
Fig. [8] 

Since the mean density profiles are similar in the 3D and ID 
cases, the lower stream densities in 3D lead to a larger number of 
streams at each radius. Counting the number of streams in ID is 
straightforward, given the similarity solution. In 3D we can esti- 
mate the number of streams crossing a given radial bin as a suit- 
able average over the stream densities of the particles it contains 
A^strcams ~ pij) {l/ps)r, where p{r) is the mean mass density in 
the bin and {l/ps)r is the average of the reciprocal of the stream 
densities of the particles. The comparison is shown in Fig. [12] The 
greatly enhanced mixing in 3D is clearly seen in this figure. At 
1% of the turnaround radius there are of order lO*" streams in the 
3D simulation but only of order 100 in the FG solution. Estimates 
of the number of streams near the Sun based on ID models (e.g. 
iNataraian & Sikiviel2005l) give severe underestimates of the expec- 
tation in more realistic models. 

In Fig. [13] we show as a function of radius the median and 
quartiles of the density at caustic crossing for all particles in the 3D 
simulation that crossed a caustic in the 0.014 Gyr immediately pre- 
ceding to - For comparison, we also show the mean density profile 
of the halo in the ID and 3D case. In the outer regions the density at 
caustic crossing is typically a thousand times the local mean halo 
density. Thus one might expect caustics to be visible in the halo 
annihilation signal at these radii. In the inner halo, however, the 
densities at caustic crossing do not rise as in the ID similarity so- 
lution (also shown for comparison) but rather stay constant or even 
drop somewhat. Below r/rta = 0.03 the median density at caustic 
crossing is smaller than the local mean halo density and one may 
expect caustics to play no significant role in the annihilation signal. 

In the inner regions, stream and caustic densities are much 
lower in the 3D simulation than in the ID similarity solution, but 
the number of caustics is only a few times greater. As a result less 
caustic-related annihilation radiation is expected in the 3D case. In 
Fig[T4]we show the ratio of intra-stream annihilation rate (i.e. the 
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Figure 15. Intra-stream contiibution to the annihilation rate as a function 
of the present-day neutralino velocity dispersion in unclustered regions and 
for different radial ranges (distinguished by colour; see legend which gives 
the range in units of the turnaround radius). The caustic contribution dom- 
inates and increases only logarithmically with decreasing dispersion. The 
dashed lines are analytic fits to the numerical results using Eq. (TJ- The 
fitting parameters a, b are shown in the brackets in the legend. 



contribution from particles which are both part of the same stream) 
to the annihilation rate estimated from the smoothed total density. 
The latter is calculated by SPH techniques using 64 neighbours. 
Figim shows that intra-stream annihilation (which is dominated 
by caustic annihilation) does not contribute strongly to the over- 
all annihilation rate. This is especially true in the inner region of 
the halo, where the stream and caustic densities are strongly sup- 
pressed. We find that for the radial interval r/rta G (0.01, 0.5) the 
intra-stream contribution is only 4%. If we focus on the outer re- 
gions we find a contribution of 24% for r/rta G (0.1, 0.5) and 
64% for r/rta G (0.2,0.5), the range where the particles turn 
around for the second time and produce the outermost caustics. 

All the above results are based on a present-day neutralino 
velocity dispersion of ah(to) = 0.03 cm/s in unclustered re- 
gions. The dominant contribution to the intra-stream annihilation 
rate comes from caustics if the velocity dispersion is small. The 
densities of caustics are proportional to l/y^Ob, and so are easily 
scaled to an arbitrary velocity dispersion. In Fig.[T5]we show the ra- 
dial dependence of the intra-stream annihilation rate for a range of 
velocity dispersions (Tb(io). These results are based on 64'^ particle 
simulations with a softening length of 0.01 rta. We plot the intra- 
stream contribution for three different radial ranges as a function 
of velocity dispersion. A good analytic fit to the numerical results 
(dashed lines) for each radial interval is given by 



(intra) 



a + b log 



1 cm/s 



(7) 



(smooth) o-b{to) 

We note that this logarithmic behaviour is expected given the anal- 
ysis of WV and the fact that the dominant contribution comes from 
caustics. From Fig. \T5\ we conclude that decreasing the velocity 
dispersion by 10 orders of magnitude from our standard value only 
increases the intra-stream annihilation contribution by a factor of 
about two. In the outer part of the halo, this would boost the total 
annihilation rate by almost a factor of two for the lowest disper- 
sion shown in the plot, but in the inner regions the caustic contri- 
butions would still be small. For a typical neutralino with a mass 



around 100 Ge V c ^ we do not confirm the substantial boost fac- 
tors claimed bv lHogai] i200lh . 



4 CONCLUSIONS 

We have use d the geodesic deviation e quation formalism recently 
developed bv lVogelsberger et al.l ( l2008h tog ether with the fully gen- 
eral tr eatment of annihilation in caustics bv lWhite & Vogelsbergerl 
( I2OO9I) to understand how fine-grained phase-space structure affects 
the annihilation radiation from the dark matter haloes that form 
from self-similar, spherical initial conditions. Such objects do not 
evolve according to the well-known one-dimensional similarity so- 
lutions. Rather, they turn into highly elongated bars as a result of 
the radial orbit instability. These bars have mean mass density and 
total velocity dispersion profiles which are very similar to those 
of the relevant similarity solutions, but the loss of spherical sym- 
metry results in orbits that fill a 3D volume and along which the 
stream density typically decreases as l/t^ rather than as 1 /t as in 
the similarity solution. At any given time, typical stream densities 
decrease slightly towards the centre of our simulation, whereas they 
increase strongly in the similarity solution. As a consequence, there 
are many more streams in the inner regions of the simulation than 
in the similarity solution. At 1% of the turnaround radius we find 
^ streams in the simulation but only about 100 in the similarity 
solution. This contradicts recent claims that th e number of streams 
near the Sun should be relatively small (e.g. iNataraian & Sikivid 
I2OO5I) but agrees with the estimates of iHelmi et al.l j2003h . The 
number of caustics changes much less dramatically between the 
two cases, with a few times more caustics in the inner regions of 
the simulation than in the comparable region of the similarity so- 
lution. This is as expected given the higher dimensionality of the 
simulation orbits. 

The impact of caustics on the annihilation signal depends on 
their density and number. Caustic densities in our simulation are 
much smaller than in the similarity solution, but their abundance 
is only modestly increased. As a result, annihilation radiation from 
caustics is less important in 3D than in ID. For example, within the 
radial range r/rta G (0.01, 0.5) caustics contribute only 4% of the 
smooth annihilation signal. If we focus on the region containing the 
outermost caustics, r/rta G (0.2, 0.5), this ratio is 64%, similar to 
that predicted by the similarity solution. Decreasing the velocity 
dispersion by 10 orders of magnitude from our standard value only 
increases the caustic contribution to the annihilation luminosity by 
a factor of about two. This is because the annihilation signal from 
caustic crossing depends only logarithmically on the dark matter 
velocity dispersion. 

Our results are based on a simplified and unrealistic halo 
formation model. However, haloes growing from ACDM initial 
conditions are expected to mix even more efficiently, because the 
fully three-dimensional character of orbits is retained and small- 
scale structure is expected to enhance the stretching of the phase- 
space sheets and hence to result in even greater dilution of their 3- 
densities. Thus, caustics will likely be less important in the ACDM 
case than in the simple isolated halo model discussed in the present 
paper. This strengthens our conclusion that ID similarity solutions 
are inadequate and misleading models for the evolution of the fine- 
grained structure of real dark matter haloes. This applies not only to 



the original FG solutio ns but also to spher i cally s ymmetric gener- 
alisations of them (e.g. ISikivie et al. l ll995Lll997l : lDuffv & Sikivid 
2008). These still give qualitatively incorrect predictions for the 
dynamical dilution of stream densities. The inclusion of baryons 
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would affect the dynamics of the dark matter in the inner halo, but 
would neither change the dimensionality of the orbits nor substan- 
tially modify their characteristic timescales. Thus no qualitative 
changes in behaviour are expected. The somewhat shorter orbital 
timescales produced by compression of the inner halo are likely, 
if anything, to accelerate mixing. Our general conclusions should 
thus be unaffected. 

Our analysis assumes that the annihilation cross-section does 
not depend on the dark matter phase-space structure. In general this 
assumption is correct, but it does not hold for a recently proposed 
mechanism which adopts an additional attractive force between 
dark matter particles in order to increase the annihilation cross- 
section t hrough the well-known Sommerfeld enhancement pro- 
cess ( e g.[Sommerfeldll93l] ; |Hisano et al.ll2004l2005l:lcirelli et all 



l2007l : lArkani-Hamed et alj 120091 ; iLattanzi & Silkl llOOsI) . Awav 
from resonances and before saturation takes place, the increase in 
the annihilation cross-section scales like 1/v, where v is the par- 
ticle encounter velocity. As a result, annihilation radiation from 
low-velocity-dispersion regions like subhaloes is enhanced as the 
inverse of the velocity dispersion (aga in as long a s this disper- 
sion exceeds the saturation level) (e.g. IBov"vI|2009|) . Conversely, 
annihilation radiation from caustics is suppressed, because the ve- 
locity dispersion there is very high. In such a scenario the radia- 
tion from individual fine-grained streams is Sommerfeld enhanced 
away from caustics because of their low internal velocity dispersion 
there. Thus the intra-stream annihilation is no longer dominated by 
caustics as it is in the standard case we have treated in the bulk of 
this paper. We will discuss this in more detail in a future publica- 
tion. 
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APPENDIX A: GDE ANALYSIS OF THE ID 
SELF-SIMILAR INFALL 

We start from the equations of motion in the similarity variables 
A = r/rta and r = t/tta introduced by EG (tta and rta are the 
turnaround time and radius of the particle under consideration) 



dr2 



V^/f^*^) . /A 
A 



(Al) 



where M is the dimensionless enclosed mass (seeEq. (0), A(r) = 
r® and 9 = 2/3 + 2/(9e). The initial conditions are A(l) = 
1 and dA/dr(l) = 0, corresponding to the Lagrangian physical 
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Figure Al. Caustic density for tlie ID similarity solution witli e = I. Red 
shows the densities using the exact d^A/dv^ values from the solution, 
while the blue line shows the result from using the Galilean-invariant esti- 
mate discussed in the text. The green line is the result obtained by MS. The 
agreement of our results with those presented in MS shows that our different 
approaches to calculating densities at caustic passage give the same result. 



coordinates q = rta (radial distance) and p — (radial velocity). 
The physical distortion tensor (see VWHS) can be related to the 
equivalent tensor expressed in terms of similarity variables through 

t 



Drp — —Dx\' 

T 



(A2) 



Daaq and Dxx[^ are two solutions to the the geodesic deviation 
equation in two-dimensional phase-space 



where D = Daao or D = D^x'^- The initial conditions are 
5(1) = = OforDAAn and 5(1) = 0,5'(1) = 1 for 

DaaJi ■ The two missing tensor components are given by r deriva- 
tives of D. To determine the stream density as a function of r we 



need the linear and second order terms 
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where A{q, p) = p — V{q) and rta is the current turnaround radius. 
The stream density for a particle is the integral of the phase-space 
density over velocity space (see WV). Using the linear and second 
order results from above, this yields 
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for the central particle of a phase-sheet, where the factor 1/A^ is 
due to the l/(siS2) prefactor associated with the stretching of the 
non-caustic directions and ^5^1/4 is the modified Bessel function. 
The intra-stream annihilation rate contribution of a small mass el- 
ement AMi is given by dPintra = Ps{t) AMi. Integrating over all 
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Figure A2. Top panel: Density at caustic crossing in units of the cosmic 
mean density and for different velocity dispersions. Density scales like 
at each caustic. The dispersions differ by factors of 10 so that the 
caustic densities differ by factors of v'TO. Bottom panel: Cumulative intra- 
stream annihilation rate counted outside to inside for the same set of dis- 
persions. Although the dispersions change by 8 orders of magnitude, the 
intra-stream annihilation rate only changes by a factor of 2 over the radial 
range shown. This is due to the logarithmic (Tt, dependence of the annihila- 
tion contribution from caustic crossing. 



mass elements yields the total intra-stream annihilation rate 
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where A/ta is the turnaround mass. 

In Fig. lAll we show the densities of the first 10 caustics calcu- 
lated using this GDE approach. For comparison to previous work 
we have chosen e = 1 and rta (to) = 1400 kpc. We used both the 
exact second-order term of the similarity solution and an estimate 
based on Galilean-invariant first-order quantities (see WV) 



(A7) 



The green dashed line shows the results of MS. All results show 
good agreement. In Fig. lA2l (top panel) we show the caustic density 
for velocity dispersions ranging from 10^ cm/s down to 10~® cm/s. 
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The maximum density scales as 1/ ,JaZ (see WV) so that all lines 
in this plot are separated by a factor of \/IO. The bottom panel 
shows the corresponding cumulative intra-stream annihilation rate 
(outside to inside). Although the dispersion varies over eight or- 
ders of magnitude, the total intra-stream contribution only changes 
by about a factor of 2. This is a consequence of the logarithmic di- 
vergence of the intra-stream annihilation luminosity near caustics 
(see WV). 



